# The Tournament of Party Strategies
# version 2.3
# james fowler
# uc san diego
# 6/7/2006
#
# vn = number of voters
# sims = number of simulations
# efreq = frequency of elections
# birth.set = which strategies (by "ptype" number) to include in tournament -- default is all
# pthresh = party minimum threshold
# pthreshperiods = number of periods a party must be below threshold to die
# starttype = starting party type -- if 0, random, otherwise the given type
# pstepsize = step size for party strategies
# plot = T/F variable to turn on/off plot of party strategies
# writeall = T/F variable to record T = all poll and election data or F = just election data
# file = name of file to write results to -- two files are recorded one with 
#        the name given and one with an "a" 
#        as a prefix which records aggregate information
# 
# party types
#
# Laver types
#  1 = STICKER
#  2 = AGGREGATOR
#  3 = HUNTER
#  4 = PREDATOR
# Entrant types
#  5 = AVERAGE
#  6 = AVOIDER
#  7 = BIGTENT
#  8 = CENTER-MASS
#  9 = FISHER
# 10 = FOLLOW-THE-LEADER
# 11 = FOOL-PROOF
# 12 = GENETY
# 13 = HALF-AGGREGATOR
# 14 = INSATIABLE-PREDATOR
# 15 = JUMPER
# 16 = KQSTRAT
# 17 = MEDIAN-VOTER-SEEKER
# 18 = MOVE-NEAR-SUCCESSFUL
# 19 = ZENO
# 20 = NICHE-HUNTER
# 21 = NICHE-PREDATOR
# 22 = STICKY-HUNTER-MEDIAN-FINDER
# 23 = OVER-UNDER
# 24 = PARASITE
# 25 = PATCHWORK
# 26 = PICK-AND-STICK
# 27 = PRAGMATIST
# 28 = RAPTOR
# 29 = SHUFFLE
#

party.tournament <- function( vn=1000, sims=1000, efreq=20, birth.set=(1:29),
                              pthresh=0.1, pthreshperiods=2, starttype=0,
                              pstepsize=0.05 , plot=T, writeall=F, file=NULL) {


 #-------------------------------------------------------------------------------------------------------------------------
 # bigtent begins here
 #-------------------------------------------------------------------------------------------------------------------------

	#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~#
	# 					BigTent Support Functions and Variables						#
	#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~#

	btGrid=NULL									# start with an empty grid
	btGridX=list()								# grid x values
	btGridY=list()								# grid y values
	btGridV=list()								# grid data values
	btData=NULL									# data to reconstruct the grid
	cip.x=NULL									# list of points on the bounding circle
	cip.y=NULL

	
	#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	# bigtents attempt to find the largest area
	#
	bigtent.adapt = function(i)
	{
		eps=0.00001						# an epsilon\
		vSmall=eps/1000					# a very small number
 		pl=NULL							# list of lines			
 		thisTent=NULL					# which tent are we working on
		gSize=8							# initial grid size
		tailLength=30					# max number of location data points to save per party
		idx=period%%tailLength			# index into a circular queue
		maxI=30						# max iterations of the optimizer
		cPts=2							# number of points around the circle to attempt to optimize about
		R=100							# circle of radius whatever

 		bigtents<-which(ptype==i)
  		others<-which(ptype!=i)

		#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

		g.distance2 = function (f,t) 		(f[1]-t[1])^2+(f[2]-t[2])^2
		g.distance = function (f,t)			sqrt(g.distance2(f,t))
		g.midpoint <- function (from, to)	(to+from)/2
		g.slope <- function (from, to)		cbind((from[,2]-to[,2])/(from[,1]-to[,1]))
		g.intercept <- function(p,m)		p[,2]-m*p[,1]
		g.lines <- function(from, to)		{ m=g.slope(from,to);cbind(m,g.intercept(from,m))}
		g.circle <- function(r) 			{curve(sqrt(r^2-x^2),-r,r,add=TRUE); curve(-sqrt(r^2-x^2),-r,r,add=TRUE)}
		g.plines <- function (from,to)
		{
			m=g.slope(from,to)
			mp=g.midpoint(from,to)
			b=g.intercept(mp,-1/m)
			cbind(-1/m,b)
		}

		m.ne <-function(a,b)	((a<b-eps) | (a>b+eps))
		m.eq <- function(a,b)	((a>b-eps) & (a<b+eps))
		m.is.singular <- function(m)	{	if(!is.qr(m)) m=qr(m); m$rank!=ncol(m$qr) }
		m.is.numbers <- function(s)
		{
			for(i in length(s))
			{
				if(is.nan(s[i])) return (FALSE)
				if(is.infinite(s[i])) return (FALSE)
				if(is.na(s[i])) return (FALSE)
			}
			return (TRUE)
		}

		bt.origin = function (plist)	rbind(mean(plist[1,]),mean(plist[2,]))
		bt.weight = function () 		0.85
		bt.gridInitialValue=function()	1/gSize^2
		bt.updateCell = function(old, new) (1-bt.weight())*old+bt.weight()*new
		bt.gridFactor = function() 		1.25
		
		bt.initcip = function()
		{
			cip.x<<-NULL
			cip.y<<-NULL
			nfp=11
			th=2*pi/nfp
			for(j in 1:nfp)
			{
				cip.x<<-c(cip.x,R*cos(th*j))
				cip.y<<-c(cip.y,R*sin(th*j))
			}					
		}


		#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
		btenergy.pcount.bbox <<- function(p)
		{
			x=p[1]; y=p[2];		  	
		 	btName=myName=pname[bigtents[thisTent]] 		
		  	for(j in 1:pn) 
  			{
  				if(j!=bigtents[thisTent])
  					{
  						if(m.ne(oploc[1,j],x) | m.ne(oploc[2,j],y))
  						{ 
							f = cbind(x,y)
							t = cbind(oploc[1,j],oploc[2,j])
							pl =rbind(pl,g.plines(f,t))						
						}
					}
				}
  						

			
			# method fails outside the bounding circle, return inf
			if((x^2+y^2)>R^2) return(vSmall)
			nlines=length(pl)/2

			if(nlines>0)
			{					
				# now compute the area within the intersection points
				g.x=btGridX[[btName]]
				g.y=btGridY[[btName]]
				g.v=btGridV[[btName]]

				# gather the values that belong to us
				if(length(g.x)>0)
				{
					# ri are the points to remove, they are not on our side
					ri=NULL
					
					for (j in 1:length(g.x))
					{
						for (k in 1:nlines)
						{
							m=pl[k,1]; 
							b=pl[k,2]; 
							c= m*x+b

							x1=g.x[j]; 
							y1=g.y[j]
							
							d = m*x1+b
							
							c1=((y>(c+eps)) & (y1<(d-eps)))
							c2=((y<(c-eps)) & (y1>(d+eps)))

							if(c1|c2)
							{
								ri=c(ri,j)
								break;
							}
						}
					}
					
					# remove points not in our intersection
					if(length(ri)>0)
					{
						g.x=g.x[-ri]
						g.y=g.y[-ri]
						g.v=g.v[-ri]
					}

					if(length(g.x>0)) 
					{
						np=length(g.x);
						a=0
						for(j in 1:np) a=a+g.v[j]
						return(1/a)
					}
				}
			}
			return(vSmall)
		}
		
		#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
		# make the grid on demand for a particular bigtent party 
		#
		makeGrid = function(btName)
		{
			gX=array(0,c(gSize,gSize))
			gY=array(0,c(gSize,gSize))
			gData=array(bt.gridInitialValue(),c(gSize,gSize))

			myData=btData[,which(btData[1,]==btName)]
			gIdx=which(btGrid[1,]==btName)
			if(length(gIdx)>0) btGrid<<-btGrid[,-gIdx]

  			min.x=min(myData[3,])
  			max.x=max(myData[3,])
  			min.y=min(myData[4,])
			max.y=max(myData[4,])
			
			# set radius of bounding circle
			R<<-g.distance(c(0,0),c(max.x,max.y))*bt.gridFactor()

			gDiameter=2*R 
			gCellSize=(gDiameter*bt.gridFactor())/(gSize-1)
			gOrigin=c(0,0)

			# calculate the x,y points (this was original 2d for graphing)
			offset=(gDiameter/2)*bt.gridFactor()
			initX=-R*bt.gridFactor()
			initY=-R*bt.gridFactor()

			for(j in 1:gSize)
			{	
				gX[j,1]=initX
				gY[gSize,j]=initY

				for(k in 2:gSize)
				{
				 	gX[j,k] = gX[j,k-1]+gCellSize
				 	gY[(gSize+1)-k,j] = gY[(gSize+2)-k,j]+gCellSize
				}
			}

			# for each datapoint in the grid
			for(j in 1:length(myData[1,]))
			{
				x=myData[3,j]		#ploc.x
				y=myData[4,j]		#ploc.y
				v=myData[5,j]		#ototal/vn
				i.x=max(min(round((x+offset)/gCellSize)+1,gSize),1)
				i.y=max(min(-round((y-offset)/gCellSize)+1,gSize),1)
				gData[i.y,i.x] = bt.updateCell(gData[i.y,i.x],v)
			}
			

			# save the results of the grid caluclation
			btGridX[[btName]]<<-gX
			btGridY[[btName]]<<-gY
			btGridV[[btName]]<<-gData
			bt.initcip()

			thisGrid=rbind(	btName,
							min.x,
							max.x,
							min.y,
							max.y,
							gDiameter,
							gCellSize,
							gOrigin[1],
							gOrigin[2]
						)
			btGrid<<-cbind(btGrid,thisGrid)
			return(thisGrid)
		}

		#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
		# bigtent begins here
		#
  		if(length(bigtents)>0 & pn > 1) 
  		{

			# eliminate old btData at index idx
			rmv=which(btData[2,]==idx)
			if(length(rmv)>0) btData<<-btData[,-rmv]

  			for(m in 1:length(bigtents))
  			{
				minV=Inf; bestX=0; bestY=0
				btName=pname[bigtents[m]]
				btData<<-cbind(btData,rbind(rep(btName,pn),
											idx,
											oploc[1,],
											oploc[2,],
											ototals/vn
											))

				myGrid=NULL
				gIndex = which(btGrid[1,]==btName)
				if(length(gIndex)>0 & period > 100)
				{
					# my grid must already be setup
					myGrid=btGrid[,gIndex[1]]
					c.x=c.y=TRUE	# rebuild on error
					if(length(oploc)>0 & length(myGrid)>0)
					{
						c.x=min(oploc[1,])< myGrid[2] | max(oploc[1,])> myGrid[3]
						c.y=min(oploc[2,])< myGrid[4] | max(oploc[2,])> myGrid[5]
					}

					# but, if some point is outside the grid, we have to rebuild the grid
					if(c.x | c.y) myGrid=makeGrid(btName)
					# otherwise, just add the points to the existing grid.
					else
					{
						gData = btGridV[[btName]]
						gDiameter=myGrid[6]
						gCellSize=myGrid[7]
						offset=(gDiameter/2)*bt.gridFactor()
						for(i in 1:pn)
						{
							x=oploc[1,i]; y=oploc[2,i]; v=ototals[i]/vn
							i.x=max(min(round((x+offset)/gCellSize)+1,gSize),1)
							i.y=max(min(-round((y-offset)/gCellSize)+1,gSize),1)
							gData[i.y,i.x] = bt.updateCell(gData[i.y,i.x],v)
						}
						btGridV[[btName]]=gData
					}
				}
				else myGrid=makeGrid(btName)

				# now optimize our location
				for(k in 1:cPts)
				{					
					thisTent=m
					x=ploc[1,bigtents[m]]
					y=ploc[2,bigtents[m]]
					th=atan(y/x)+(2*pi/cPts)*(k-1)
					r=sqrt(x^2+y^2)*(1+runif(1,0,.25)*(-1)^(period%%2)*(k>1))#allows testing random points in an annular region, based on own radius
					x=r*cos(th)
					y=r*sin(th)

					#adapt					
					if(pn>1)
					{
						try(
						{
							result=optim(c(x,y),
									try(btenergy.pcount.bbox,silent=TRUE),
									method="BFGS",
									control=list(reltol=0.0001,maxit=maxI)
									)

							if(result$value < minV)
							{
								bestX=result$par[1]
								bestY=result$par[2]						
								minV=result$value
							}
						},silent=TRUE)
					}
				}

				if(is.finite(minV))
				{
   					ploc[1,bigtents[m]]<<-bestX
   					ploc[2,bigtents[m]]<<-bestY		
				}
			}
  		}  	
	}

 #-------------------------------------------------------------------------------------------------------------------------
 # bigtent ends here
 #-------------------------------------------------------------------------------------------------------------------------


 # initialize strategy list and associated colors

 plabels<-c("STICKER","AGGREGATOR","HUNTER","PREDATOR",
            "AVERAGE","AVOIDER","BIGTENT","CENTER-MASS","FISHER",
            "FOLLOW-THE-LEADER","FOOL-PROOF","GENETY","HALF-AGGREGATOR",
            "INSATIABLE-PREDATOR","JUMPER","KQSTRAT","MEDIAN-VOTER-SEEKER",
            "MOVE-NEAR-SUCCESSFUL","ZENO","NICHE-HUNTER",
            "NICHE-PREDATOR","STICKY-HUNTER-MEDIAN-FINDER","OVER-UNDER","PARASITE",
            "PATCHWORK","PICK-AND-STICK","PRAGMATIST","RAPTOR","SHUFFLE")
 pcolors<-rainbow(29,gamma=0.5)

 # initialize voter ideal points, poll/election indicator
 videal<-rbind(rnorm(vn),rnorm(vn))
 pestr<-c("poll","election")
 epsilon<-1/10^10
 if(plot) par(mfrow=c(1,3),las=3, omi=c(0.5,0,0,0), mai=c(1,0.4,1.5,0.1))
 
 # initialize party location, name, and type vectors
 ototals<-angle<-oploc<-pbt<-ploc<-pmloc<-pmedloc<-pname<-ptype<-ftldim<-NULL
 umlocch<-c(0,0)

 # initialize proportional vote counts
 pcount<-mcount<-rep(1/29,29)
 etop<-1
 
 # intialize datafiles if there is one
 if(length(file)>0) write.table(t(c("Election","PartyID","Type","VoteRecd",
                                "X","Y","Angle")), file=file, 
                                row.names=F,col.names=F,quote=F,sep=",")
 # initialize fisher vars
 farray=NULL
 ftoppos=NULL
 fstepsizes=NULL

 # initialize vars for kqstrat

 KQstrat14 <- function(ploc, locat1.history, locat2.history, votes.history,
                     pthresh, strat.indic){
  
  n <- ncol(locat1.history) # number of parties
  stepsize <- 1.46e-4
  lastmove.dist.vec <- sqrt((locat1.history[1,] - locat1.history[2,])^2 +
                            (locat2.history[1,] - locat2.history[2,])^2)
  lastmove.dist.vec[is.na(lastmove.dist.vec)] <- stepsize
  for (i in 1:n){
    if (strat.indic[i]){      
      if (sum(!is.na(votes.history[1:3,i])) == 3 &
          sum(votes.history[1:3,i] < pthresh) == 3){
        if (sum( (locat1.history[1,i] -
                  locat1.history[1:3,i]) < stepsize,
                na.rm=TRUE) == 3 &
            sum((locat2.history[1,i] -
                 locat2.history[1:3,i]) < stepsize,
                na.rm=TRUE) == 3
            ){
          possib.ind <- which(votes.history[1,]>pthresh &
                              lastmove.dist.vec!=stepsize)
          
          if (length(possib.ind) == 0){
            possib.ind <- which(votes.history[1,]>pthresh)
          }
          ind <- sample(possib.ind, 1)
          ploc[1,i] <- locat1.history[1,ind] + rnorm(1, m=0, s=.0625)
          ploc[2,i] <- locat2.history[1,ind] + rnorm(1, m=0, s=.0625)
        }
        else{
          a <- runif(1, 0, 2*pi)
          ploc[1,i] <- ploc[1,i] + stepsize * cos(a)
          ploc[2,i] <- ploc[2,i] + stepsize * sin(a)
        }        
      }
    }    
  }
  return(ploc)  
}


  locat1.history <- NULL
  locat2.history <- NULL
  votes.history <- NULL
  history.length <- 100 


  for(period in 1:sims) {
 
  #
  # party birth
  #
 
  if(period%%efreq==1) {
 
   # party name, failure counter, previous total counter, and angle
   pname<-c(pname,1+(period-1)/efreq)
   pbt<-c(pbt,0)
   ototals<-c(ototals,0)
   angle<-c(angle,runif(1,0,2*pi))
 
   # draw party type
   if(period==1&starttype!=0) {
    ptype<-starttype
   } else {
    ptype<-c(ptype,sample(c(birth.set,birth.set),1))
   }
 
   # count parties
   pn<-length(ptype)
 
   # party initial position
   ploc<-cbind(ploc,rbind(rnorm(1),rnorm(1)))
   pmloc<-cbind(pmloc,ploc[,pn])
   pmedloc<-cbind(pmedloc,ploc[,pn])
   oploc<-cbind(oploc,ploc[,pn])
   ftldim<-c(ftldim,sample(c(T,F),1)) #follow the leader variable

   # location histories
   locat1.history <- cbind(locat1.history, rep(NA, history.length))
   locat2.history <- cbind(locat2.history, rep(NA, history.length))      
   votes.history <- cbind(votes.history, rep(NA, history.length))
 
  }

  # update histories
  locat1.history <- rbind(ploc[1,], locat1.history)
  locat1.history <- as.matrix(locat1.history[1:history.length,])      
  locat2.history <- rbind(ploc[2,], locat2.history)
  locat2.history <- as.matrix(locat2.history[1:history.length,])      

  #
  # take poll / hold election
  #
 
  # find distances between parties and voters / respondents
  dist<-matrix(0,nrow=pn,ncol=vn)
  for(i in 1:pn) {
    
   dist[i,]<-colSums((videal-ploc[,i])^2)+ # note sqrt not taken to save operation
             runif(vn,-epsilon,epsilon)    # and epsilon added to break ties
  }
 
  # voters choose closest party (epsilon added to break ties)
  vchoice<-apply(dist,2,which.min)
  
  # election / poll results
  ptotals<-hist(vchoice,breaks=0:pn,plot=F)$counts
  votes.history <- rbind(ptotals/vn, votes.history)
  votes.history <- as.matrix(votes.history[1:history.length,])
 
  ddm<-c(median(ploc[1,vchoice]),median(ploc[2,vchoice]))

  if(period%%efreq==0) {
 
   # record counts and then randomly 
   # choose party type for next birth
   for(i in 1:29) {
    div<-sum(ptype==i)
    pcount[i]<-div
    if(div==0) div<-1
    mcount[i]<-sum(ptotals[ptype==i])/vn/div
   }
  }

  #
  # display information
  #
 
  # print party info
  if(writeall||period%%efreq==0) {
   if(length(file)==0) {
    write.table(paste("Period",period,pestr[(period%%efreq==0)+1]),
     row.names=F,col.names=F,quote=F)
    write.table(format(rbind(c("Election","Party ID","Type","Vote Recd",
     "X","Y","Angle"),
     cbind(rep(period/efreq,pn),pname,ptype,round(ptotals/vn,3),
     round(t(ploc),3),round(angle,3)))),
     row.names=F,col.names=F,quote=F)
   } else {
    write.table(cbind(rep(period/efreq,pn), pname, ptype, round(ptotals/vn,3),
     round(t(ploc),3), round(angle,3)), file=file, row.names=F, col.names=F, 
     append=T, sep=",")
    write.table(t(c(mcount,pcount)), file=paste("a",file,sep=""),
     row.names=F, col.names=F, append=T, sep=",")
   }
  }

  # plot parties
  if(plot) {
   plot(0,0,type="n",xlim=c(-3,3),ylim=c(-3,3),xlab="",ylab="",main="Location")
   text(t(ploc),labels=as.character(pname),col=pcolors[ptype])
   barplot(mcount,names.arg=plabels,cex.names=0.8,
    col=pcolors,main="Votes per Party")
   barplot(mcount*pcount,names.arg=plabels,cex.names=0.8,
    col=pcolors,main="Total Votes")
  }  

  #
  # calculate available information
  #
 
  # mean and median party members' location
  for(i in 1:pn) {
   pslist<-which(vchoice==i)
   if(length(pslist)==0) {
    pmedloc[,i]<-pmloc[,i]<-ploc[,i]
   }
   if(length(pslist)==1) {
    pmedloc[,i]<-pmloc[,i]<-videal[,pslist]
   }
   if(length(pslist)>1) {
    pmloc[,i]<-rowMeans(videal[,pslist])
    pmedloc[1,i]<-median(videal[1,pslist])
    pmedloc[2,i]<-median(videal[2,pslist])
   }
  }
 
  # mean of all party locations
  if(pn==1) {
   mloc<-ploc 
  } else {
   mloc<-colSums(ptotals*t(ploc))/vn  
  }

  #
  # party death
  #
 
  if(period%%efreq==0) { 
 
   # which parties fell below threshold?
   pweak<-as.numeric(ptotals<pthresh*vn)
 
   # update number of times party falls below threshold
   pbt<-(pbt+pweak)*pweak
 
   # kill parties that have been below threshold for too long
   pdead<-which(pbt==pthreshperiods)
   if(length(pdead)>0) {
    pn<-pn-length(pdead)
    pname<-pname[-pdead]
    ptype<-ptype[-pdead]
    ploc<-ploc[,-pdead]
    pmloc<-pmloc[,-pdead]
    pmedloc<-pmedloc[,-pdead]
    oploc<-oploc[,-pdead]
    angle<-angle[-pdead]
    pbt<-pbt[-pdead]
    ptotals<-ptotals[-pdead]
    ototals<-ototals[-pdead]
    ftldim<-ftldim[-pdead]
    locat1.history <- locat1.history[,-pdead]
    locat2.history <- locat2.history[,-pdead]
    votes.history  <- votes.history[,-pdead]

    #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	#   bigtent dies
	#
	# eliminate dead btData
	for(j in 1:length(pdead))  btData<<-btData[,-which(btData[1,]==pname[j])]
    
	# end of bigtent death
    #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

   }
 
  }
 
  
  #
  # calculate more available information
  #
 
  # angle and distance to top party
  prank<-rank(ptotals,ties.method="random")
  top<-which(prank==pn)
  topploc<-ploc[,top]
  angletop<-atan((ploc[2,top]-ploc[2,])/(ploc[1,top]-ploc[1,]))+
            pi*(ploc[1,top]<ploc[1,])%%(2*pi)
  angletop[which(is.nan(angletop))]<-runif(sum(is.nan(angletop)),0,2*pi)
  disttop<-sqrt(colSums((ploc[,top]-ploc)^2))
  steptop<-apply(matrix(c(disttop,rep(pstepsize,pn)),ncol=2),1,min)

  if(period%%efreq==0) etop<-top
  etopploc<-ploc[,etop]

  # location of 2nd party (or 1st if second doesn't exist)
  second<-which(prank==max(pn-1,1))
  anglesecond<-atan((ploc[2,second]-ploc[2,])/(ploc[1,second]-ploc[1,]))+
            pi*(ploc[1,second]<ploc[1,])%%(2*pi)
  anglesecond[which(is.nan(anglesecond))]<-runif(sum(is.nan(anglesecond)),0,2*pi)

  # unique location of 2nd party
  usecond<-top
  if(pn<3|ploc[1,top]!=ploc[1,second]|ploc[2,top]!=ploc[2,second]) {
   usecond<-second
  } else {
   ius<-pn-1
   while(ius>1) {
    ius<-ius-1
    utest<-which(prank==ius)
    if(ploc[1,top]!=ploc[1,utest]|ploc[2,top]!=ploc[2,utest]) { 
     usecond<-utest
     ius<-1
    }
   }
  }

  # angle and distance to weighted mean of parties 
  anglem<-atan((mloc[2]-ploc[2,])/(mloc[1]-ploc[1,]))+
            pi*(mloc[1]<ploc[1,])%%(2*pi)
  anglem[which(is.nan(anglem))]<-runif(sum(is.nan(anglem)),0,2*pi)
  distm<-sqrt(colSums((mloc-ploc)^2))
 
  # angle and distance to ddm
  angled<-atan((ddm[2]-ploc[2,])/(ddm[1]-ploc[1,]))+
            pi*(ddm[1]<ploc[1,])%%(2*pi)
  angled[which(is.nan(angled))]<-runif(sum(is.nan(angled)),0,2*pi)

  # store this round's recent ploc
  o2ploc<-ploc

 
  #
  # parties adapt
  #
 
  # 1.stickers stick
 
  # 2.aggregators move to party mean
  ploc[,ptype==2]<-pmloc[,ptype==2]
  
  # 3.hunters win or tie - stay the course, lose - shift direction
  hunters<-which(ptype==3)
  if(length(hunters)>0) {
   angle[hunters]<-(angle[hunters]+(ototals[hunters]>=ptotals[hunters])*
                    runif(length(hunters),pi/2,3*pi/2))%%(2*pi)
   ploc[1,hunters]<-ploc[1,hunters]+pstepsize*cos(angle[hunters])
   ploc[2,hunters]<-ploc[2,hunters]+pstepsize*sin(angle[hunters])
  }
 
  # 4.predators move to largest party (unless they are the largest)
  predators<-which(ptype==4&prank!=pn)
  if(length(predators)>0) {
   ploc[1,predators]<-ploc[1,predators]+
    steptop[predators]*cos(angletop[predators])
   ploc[2,predators]<-ploc[2,predators]+
    steptop[predators]*sin(angletop[predators])
  }
 
  # 5.average move to average party location weighted by votes
  averages<-which(ptype==5)
  nn<-length(averages)
  if(nn>0) {
   ploc[,averages]<-matrix(rep(mloc,nn),nrow=2)
  }
 
  # 6.avoider move to outside average party location weighted by votes
  avoiders<-which(ptype==6)
  nn<-length(avoiders)
  if(nn>0) {
   aa<-runif(nn,-1,1)
   bb<-sample(c(-1,1),nn,replace=T)
   cc<-sample(c(-1,1),nn,replace=T)
   ploc[,avoiders]<-matrix(rep(mloc,rep(nn,2)),ncol=2)+
                    cbind(aa*bb,sqrt(1-aa)*cc)
  }

  # 7. bigtent adapts

  try(bigtent.adapt(7))

  # 8.centermass moves 1/20 towards average party location weighted by votes
  centermasses<-which(ptype==8)
  if(length(centermasses)>0) {
   ploc[1,centermasses]<-ploc[1,centermasses]+distm[centermasses]*cos(anglem[centermasses])/20
   ploc[2,centermasses]<-ploc[2,centermasses]+distm[centermasses]*sin(anglem[centermasses])/20
  }

  # 9.fishers do random search and then choose good spot before election
  fishers<-which(ptype==9)
  numfishers=length(fishers)

  if(numfishers>0) {

    # Initialize variables during the period after each election
    if (period%%efreq==1) {
      # fstepsizes contains the size of the steps a party takes when fishing
      # - varies for each party according to past success
      fstepsizes=NULL

      # farray stores the each fisher party's x and y position, and level of
      # voter support for each period during an election cycle of efreq periods
      farray=array(0, c(efreq,numfishers,3))
    }

    # store polling data and position for each fisher party for every period
    # between elections
    for (fisherid in 1:numfishers) {
      farray[period%%efreq,fisherid,1]=ploc[1,fishers[fisherid]]
      farray[period%%efreq,fisherid,2]=ploc[2,fishers[fisherid]]
      farray[period%%efreq,fisherid,3]=ptotals[fishers[fisherid]]
    }
 
    if (period%%efreq==0) {
      # don't move, just recapture voter support for changed environment
      # (new party entry).  In future may just use old total and test a
      # new position here

    } else if (period%%efreq==1) {
      for (fisherid in 1:numfishers) {

        if (ototals[fishers[fisherid]]>=pthresh*vn) {
          # if the party's vote share in the last election was above the
          # threshold for party survival, use a smaller step size during
          # polling before the next election and start polling in new positions
          fstepsizes=c(fstepsizes, 0.2)
          angle[fishers[fisherid]]<-(angle[fishers[fisherid]]+(1)*runif(1,pi/4,3*pi/4))%%(2*pi)
          ploc[1,fishers[fisherid]]<-ploc[1,fishers[fisherid]]+fstepsizes[fisherid]*cos(angle[fishers[fisherid]])
          ploc[2,fishers[fisherid]]<-ploc[2,fishers[fisherid]]+fstepsizes[fisherid]*sin(angle[fishers[fisherid]])
        } else {
          # otherwise use a larger step size, then start polling from the approximated
          # mean since the party could be born far from voters and so that the party
          # does not stray too far from the mean initially.
          fstepsizes=c(fstepsizes, 1)
         
          # approximate the mean of all voters to position new or struggling fishers
          # near the mean voter position
          xmean=0
          ymean=0
          for (pid in 1:pn) {
            xmean=sum(xmean,ptotals[pid]*ploc[1,pid])
            ymean=sum(ymean,ptotals[pid]*ploc[2,pid])
          }
          xmean=xmean/vn
          ymean=ymean/vn
          ploc[1,fishers[fisherid]]=xmean
          ploc[2,fishers[fisherid]]=ymean
        }
      }
    } else if (period%%efreq<(3/5)*efreq) {

      # Continue to use the first 3/5 of the period to learn about
      # voter preferences by *randomly* casting around the voter environment
      # for policy positions and storing voter support levels (stored by code above)
      for (fisherid in 1:numfishers) {
        angle[fishers[fisherid]]<-(angle[fishers[fisherid]]+(1)*runif(1,pi/4,3*pi/4))%%(2*pi)
        ploc[1,fishers[fisherid]]<-ploc[1,fishers[fisherid]]+fstepsizes[fisherid]*cos(angle[fishers[fisherid]])
        ploc[2,fishers[fisherid]]<-ploc[2,fishers[fisherid]]+fstepsizes[fisherid]*sin(angle[fishers[fisherid]])
      }

    } else if (period%%efreq<(4/5)*efreq) {
      # Use the next 1/5 of the period to revist the top positions found in
      # the random search and re-rank those positions according to the latest
      # polling.

      # Rank the top positions
      if (period%%efreq==(3/5)*efreq) {
        ftoppos=array(0, c(efreq/5,numfishers,2)) # x and y loc
        for (fisherid in 1:numfishers) {
          for (rank in 1:(efreq/5)) {
            topsupport=0
            topsequence=0
            for (sequence in 1:((3/5)*efreq)) {
              if ((farray[sequence,fisherid,3]>topsupport)) {
                ftoppos[rank,fisherid,1]=farray[sequence,fisherid,1]
                ftoppos[rank,fisherid,2]=farray[sequence,fisherid,2]
                topsupport=farray[sequence,fisherid,3]
                # used to zero out when all is done
                topsequence=sequence
              }
            }
            # zero out the top supported position for this fisher
            # so it is not reused
            farray[topsequence,fisherid,3]=0
          }
        }
      } 
      
      # Then revist these top positions
      # farray will store their voter support totals in earlier code
      for (fisherid in 1:numfishers) {
        ploc[1,fishers[fisherid]]=ftoppos[(period%%efreq)-((3/5)*efreq)+1,fisherid,1]
        ploc[2,fishers[fisherid]]=ftoppos[(period%%efreq)-((3/5)*efreq)+1,fisherid,2]
      }
      
    } else if (period%%efreq==(4/5)*efreq) {

      # Rank the polling totals gathered in farray from the previous moves
      ftoppos=array(0, c(efreq/5,numfishers,2)) # x and y loc
      for (fisherid in 1:numfishers) {
        for (rank in 1:(efreq/5)) {
          topsupport=0
          topsequence=0
          for (sequence in ((3/5)*efreq+1):((4/5)*efreq)) {
            if ((farray[sequence,fisherid,3]>topsupport)) {
              ftoppos[rank,fisherid,1]=farray[sequence,fisherid,1]
              ftoppos[rank,fisherid,2]=farray[sequence,fisherid,2]
              topsupport=farray[sequence,fisherid,3]
              # used to zero out when all is done
              topsequence=sequence
            }
          }
          # zero out the top supported position for this fisher
          # so it is not reused
          farray[topsequence,fisherid,3]=0
        }
        # set the fisher to the best position determined above
        ploc[1,fishers[fisherid]]=ftoppos[1,fisherid,1]
        ploc[2,fishers[fisherid]]=ftoppos[1,fisherid,2]
      }
   
    } else if (period%%efreq<efreq-1) {
      for (fisherid in 1:numfishers) {
        if ((period%%efreq==(4/5)*efreq+1) && (ptotals[fishers[fisherid]]<pthresh*vn)) {
          # if the vote share is less than the threshhold for parties
          # return to the second highest ranked position
          ploc[1,fishers[fisherid]]=ftoppos[2,fisherid,1]
          ploc[2,fishers[fisherid]]=ftoppos[2,fisherid,2]
        } else {
          # "fish" until the just before the election
          angle[fishers[fisherid]]<-(angle[fishers[fisherid]]+(1)*runif(1,pi/4,3*pi/4))%%(2*pi)
          ploc[1,fishers[fisherid]]<-ploc[1,fishers[fisherid]]+(0.05)*cos(angle[fishers[fisherid]])
          ploc[2,fishers[fisherid]]<-ploc[2,fishers[fisherid]]+(0.05)*sin(angle[fishers[fisherid]])
        }
      }

    } else {
      # for the election go to the highest ranked position in the last 1/5 efreq - 1
      for (fisherid in 1:numfishers) {
        topsupport=0
        fxloc=0
        fyloc=0
        # Find the recently-visited position with the most voter support
        for (sequence in ((4/5)*efreq):(efreq-2)) {
          if ((farray[sequence,fisherid,3]>topsupport)) {
            fxloc=farray[sequence,fisherid,1]
            fyloc=farray[sequence,fisherid,2]
            topsupport=farray[sequence,fisherid,3]
          }
        }
        # return to the position that earned the most support
        ploc[1,fishers[fisherid]]=fxloc
        ploc[2,fishers[fisherid]]=fyloc
      }
    }

  }


  # 10.followtheleaders move between two top vote getters and then move toward top on one dim.
  followtheleaders<-which(ptype==10&prank!=pn)
  if(period%%efreq==1|period%%efreq==2) followtheleaders<-followtheleaders[followtheleaders!=pn]
  if(length(followtheleaders)>0) {
   ftli<-intersect(which(ptotals<=ototals),followtheleaders)
   ftldim[ftli]<-!ftldim[ftli]
   ploc[1,followtheleaders]<-ploc[1,followtheleaders]+
                             0.05*sign(cos(angletop[followtheleaders]))*
                             (!ftldim[followtheleaders])
   ploc[2,followtheleaders]<-ploc[2,followtheleaders]+
                             0.05*sign(sin(angletop[followtheleaders]))*
                             ftldim[followtheleaders]
  }
  if((ptype[pn]==10)&(period%%efreq==1)&top!=pn) ploc[,pn]<-(ploc[,top]+ploc[,second])/2
  if((ptype[pn]==10)&(period%%efreq==2)&top!=pn) {
   ftldim[pn]<-abs(ploc[2,pn]-ploc[2,top])>abs(ploc[1,pn]-ploc[1,top])
   ploc[1,pn]<-ploc[1,pn]+0.05*sign(cos(angletop[pn]))*(!ftldim[pn])
   ploc[2,pn]<-ploc[2,pn]+0.05*sign(sin(angletop[pn]))*ftldim[pn]
  }

  # 11. foolproof - shift in same direction as all other parties and at 1/5 speed
  foolproof<-which(ptype==11)
  if(period%%efreq==1) {
   foolproof<-foolproof[foolproof!=pn]
   if(ptype[pn]==11) ploc[,pn]<-mloc
  }
  if(length(foolproof)>0&pn>1) {
   ploc[1,foolproof]<-ploc[1,foolproof]+0.2*(umlocch[1]*pn
                      -(ploc[1,foolproof]-oploc[1,foolproof]))/(pn-1)
   ploc[2,foolproof]<-ploc[2,foolproof]+0.2*(umlocch[2]*pn
                      -(ploc[2,foolproof]-oploc[2,foolproof]))/(pn-1)
  }

  # 12. Genety - use a genetic algorithm
  genety <- which(ptype==12)
  genlen <- length(genety)
  w2     <- c(0.200939974226058, 1.37194332918152,  0.18451363844797);

  if(genlen>0) {

    orest <- rowSums(ploc*rbind(ptotals,ptotals))/vn;
    dorig <- -sqrt((ploc[1,genety]-orest[1])^2 + (ploc[2,genety]-orest[2])^2);
    aorig <- atan((ploc[2,genety]-orest[2])/(ploc[1,genety])-orest[1]) + pi*(ploc[1,genety]<orest[1])%%(2*pi);

    dmean <- -sqrt((ploc[1,genety]-pmloc[1,genety])^2 + (ploc[2,genety]-pmloc[2,genety])^2);
    amean <- atan((ploc[2,genety]-pmloc[2,genety])/(ploc[1,genety]-pmloc[1,genety])) + 
             pi*(ploc[1,genety]<pmloc[1,genety])%%(2*pi)

    dhunt <- rep(0.5,genlen);
    ahunt <- (angle[genety] + (ototals[genety]>=ptotals[genety])*runif(1,pi/2,3*pi/2))%%(2*pi)

    ds <- rbind(dorig, dmean, dhunt);
    as <- rbind(aorig, amean, ahunt);
    as[is.nan(as)] <- 0

    temploc <- cbind(ploc[,genety]);
    w <- t(t(matrix(1,3,genlen)*w2));
    ploc[1,genety] <- ploc[1,genety] + colSums( w * cbind(ds) * cbind(cos(as)) );
    ploc[2,genety] <- ploc[2,genety] + colSums( w * cbind(ds) * cbind(sin(as)) );

    angle[genety]  <- atan((ploc[2,genety]-temploc[2,])/(ploc[1,genety]-temploc[1,])) + 
                      pi*(ploc[1,genety]<temploc[1,])%%(2*pi);
    angle[is.nan(angle[genety])]  <- 0;
  }

  # 13. half-aggregator - 1st move aggregator - after do hunter on one fixed dimension
  halfags<-which(ptype==13)
  if(period%%efreq==1) {
   halfags<-halfags[halfags!=pn]
   if(ptype[pn]==13) {
    ploc[,pn]<-pmloc[,pn]
    angle[pn]<-sample(pi*(0:3)/2,1)
   }
  }
  if(period%%efreq==2) halfags<-halfags[halfags!=pn]
  if(length(halfags)>0) {
   angle[halfags]<-(angle[halfags]+
                   (ototals[halfags]>=ptotals[halfags])*pi)%%(2*pi)
   ploc[1,halfags]<-ploc[1,halfags]+0.05*cos(angle[halfags])
   ploc[2,halfags]<-ploc[2,halfags]+0.05*sin(angle[halfags])
  }
  if((ptype[pn]==13)&(period%%efreq==2)) {
   ploc[1,halfags]<-ploc[1,halfags]+0.05*cos(angle[halfags])
   ploc[2,halfags]<-ploc[2,halfags]+0.05*sin(angle[halfags])
  }

  # 14. insatiable predators move to largest party (unless they are the largest - move to second)
  ipredators<-which(ptype==14&prank!=pn)
  if(length(ipredators)>0) {
   ploc[1,ipredators]<-ploc[1,ipredators]+
    steptop[ipredators]*cos(angletop[ipredators])
   ploc[2,ipredators]<-ploc[2,ipredators]+
    steptop[ipredators]*sin(angletop[ipredators])
  }
  if(ptype[top]==14&pn>1) {
   ploc[1,top]<-ploc[1,top]+
    0.05*cos(anglesecond[top])
   ploc[2,top]<-ploc[2,top]+
    0.05*sin(anglesecond[top])
  }
 
  # 15. jumper move opposite of first party on other side of second
  jumpers<-which(ptype==15)
  if(length(jumpers)>0) {
    ploc[1,jumpers]<-ploc[1,top]+((ploc[1,top]>ploc[1,usecond])*2-1)*disttop[usecond]/2
    ploc[2,jumpers]<-ploc[2,top]+((ploc[2,top]>ploc[2,usecond])*2-1)*disttop[usecond]/2
  }

  # 16. KQstrat
  KQstrats <- ptype == 16
  ploc <- KQstrat14(ploc, locat1.history, locat2.history, votes.history,
                    pthresh, KQstrats)

  # 17.median finders move to ddm
  mfinders<-which(ptype==17&(ploc[1,]!=ddm[1]|ploc[2,]!=ddm[2]))
  if(length(mfinders)>0) {
   ploc[1,mfinders]<-ploc[1,mfinders]+
    0.05*cos(angled[mfinders])
   ploc[2,mfinders]<-ploc[2,mfinders]+
    0.05*sin(angled[mfinders])
  }

  # 18. move near successful party
  mnsp<-which(ptype==18)
  mnspl<-length(mnsp)
  if(mnspl>0) {
   angle[mnsp]<-runif(mnspl,0,2*pi)
   if(sum(pbt==0)>0) {
    selparties<-sample(which(pbt==0),mnspl,replace=T)
    ploc[1,mnsp]<-o2ploc[1,selparties]+0.3*cos(angle[mnsp])
    ploc[2,mnsp]<-o2ploc[2,selparties]+0.3*sin(angle[mnsp])
   }
  }

  # 19.zenos move half distance to largest party (unless they are the largest)
  zenos<-which(ptype==19&prank!=pn)
  if(length(zenos)>0) {
   ploc[1,zenos]<-ploc[1,zenos]+
    0.5*disttop[zenos]*cos(angletop[zenos])
   ploc[2,zenos]<-ploc[2,zenos]+
    0.5*disttop[zenos]*sin(angletop[zenos])
  }
 
  # 20. Niche-Hunter win or tie - stay the course, lose - move towards adjacent
  nhunters<-which(ptype==20)
  lnhunters<-length(nhunters)
  if(lnhunters>0&pn>1) {
   nhdist<-rep(0.05,lnhunters)
   for(iln in 1:lnhunters) {
    if(ototals[nhunters[iln]]>=ptotals[nhunters[iln]]) {
     anglenh<-atan((-ploc[2,nhunters[iln]]+ploc[2,])/(-ploc[1,nhunters[iln]]+ploc[1,]))+
              pi*(ploc[1,nhunters[iln]]>ploc[1,])%%(2*pi)
     anglenh[which(is.nan(anglenh))]<-runif(sum(is.nan(anglenh)),0,2*pi)
     quadrant<-trunc(anglenh/(pi/2))
     distnh<-sqrt(colSums((ploc[,nhunters[iln]]-ploc)^2))
     distnh[nhunters[iln]]<-Inf
     adjacents<-NULL
     for(jln in 0:3) {
      adjacents<-c(adjacents,(1:pn)[quadrant==jln][which.min(distnh[quadrant==jln])])
     }
     nhtarget<-adjacents[which.max(ptotals[adjacents])]
     angle[nhunters[iln]]<-anglenh[nhtarget]
     nhdist[iln]<-min(0.05,0.5*sqrt(sum((ploc[,nhunters[iln]]-ploc[,nhtarget])^2)))
    }
   }
   ploc[1,nhunters]<-ploc[1,nhunters]+nhdist*cos(angle[nhunters])
   ploc[2,nhunters]<-ploc[2,nhunters]+nhdist*sin(angle[nhunters])
  }

  # 21. Niche-Predator move to largest adjacent party - if largest stay still
  npredators<-which(ptype==21&prank!=pn)
  lnpredators<-length(npredators)
  if(lnpredators>0) {
   nhdist<-rep(0.05,lnpredators)
   for(iln in 1:lnpredators) {
    anglenh<-atan((-ploc[2,npredators[iln]]+ploc[2,])/(-ploc[1,npredators[iln]]+ploc[1,]))+
             pi*(ploc[1,npredators[iln]]>ploc[1,])%%(2*pi)
    anglenh[which(is.nan(anglenh))]<-runif(sum(is.nan(anglenh)),0,2*pi)
    quadrant<-trunc(anglenh/(pi/2))
    distnh<-sqrt(colSums((ploc[,npredators[iln]]-ploc)^2))
    distnh[npredators[iln]]<-Inf
    adjacents<-NULL
    for(jln in 0:3) {
     adjacents<-c(adjacents,(1:pn)[quadrant==jln][which.min(distnh[quadrant==jln])])
    }
    nhtarget<-adjacents[which.max(ptotals[adjacents])]
    angle[npredators[iln]]<-anglenh[nhtarget]
    nhdist[iln]<-min(0.05,0.5*sqrt(sum((ploc[,npredators[iln]]-ploc[,nhtarget])^2)))
   }
   ploc[1,npredators]<-ploc[1,npredators]+nhdist*cos(angle[npredators])
   ploc[2,npredators]<-ploc[2,npredators]+nhdist*sin(angle[npredators])
  }

  # 22. sticky-hunter/median-finder
  shmfs<-which(ptype==22)
  if(length(shmfs)>0) {
   shmfind<-which((ptotals[shmfs]<vn*0.1)&(ototals[shmfs]<vn*0.1)&pbt[shmfs]==0)
   lshmfind<-length(shmfind)
   if(lshmfind>0) {
    angle[shmfs[shmfind]]<-(angle[shmfs[shmfind]]+
                     (ototals[shmfs[shmfind]]>=ptotals[shmfs[shmfind]])*
                     (runif(lshmfind,pi/2,5*pi/6)+sample(c(0,2*pi/3),lshmfind,replace=T)))%%(2*pi)
    ploc[1,shmfs[shmfind]]<-ploc[1,shmfs[shmfind]]+0.135*cos(angle[shmfs[shmfind]])
    ploc[2,shmfs[shmfind]]<-ploc[2,shmfs[shmfind]]+0.135*sin(angle[shmfs[shmfind]])
   }
   shmfind<-which(pbt[shmfs]>0)
   if(length(shmfind)>0) {
    ploc[,shmfs[shmfind]]<-pmedloc[,shmfs[shmfind]]
   }
  }

  # 23. over-under
  overunders<-which(ptype==23)
  loverunders<-length(overunders)
  if(loverunders>0&pn>1) {
   for(iou in 1:loverunders) {
    if(period%%efreq!=19) {
     randparty<-sample(c((1:pn)[-overunders[iou]],(1:pn)[-overunders[iou]]),1)
     ploc[,overunders[iou]]<-o2ploc[,randparty]+c(20,20)
    } else {
     ploc[,overunders[iou]]<-etopploc+runif(2,0.001,0.05)*
                             (2*(angle[overunders[iou]]>pi)-1)
    }
    if(period%%efreq==19&pbt[overunders[iou]]>0) {
     angle[overunders[iou]]<-(angle[overunders[iou]]+pi)%%(2*pi)
    }
   }
  }

  # 24. parasite
  if(period%%efreq==19) {
   parasites<-which(ptype==24)
   angle[parasites]<-runif(length(parasites),0,2*pi)
   ploc[1,parasites]<-topploc[1]+0.01*cos(angle[parasites])
   ploc[2,parasites]<-topploc[2]+0.01*sin(angle[parasites])
  }

  # 25. patchwork
  patchworks<-which(ptype==25&prank!=pn)
  if(length(patchworks)>0) {
   ploc[1,patchworks]<-ploc[1,patchworks]+
    steptop[patchworks]*cos(angletop[patchworks])
   ploc[2,patchworks]<-ploc[2,patchworks]+
    steptop[patchworks]*sin(angletop[patchworks])
  }
  patchworks<-which(ptype==25&prank==pn)
  if(length(patchworks)>0) ploc[,patchworks]<-pmloc[,patchworks]

  
  # 26.Pick&Stick
  pickstick<-which(ptype==26)
  if(length(pickstick)>0) {  
   for(i in 1:pn) {
     if (period%%efreq==1) {
       indicator<-toptotals<-bestxcoord<-bestycoord<-c(1:pn)
       toptotals[i]<-0
       bestxcoord[i]<-0
       bestycoord[i]<-0 
       indicator[i]<-0
     }
     if ((pname[i]-1)*20+19 >= period){
       if (ptype[i]==26){
         if (ptotals[i]>toptotals[i]) {
           toptotals[i]<-ptotals[i]
           bestxcoord[i]<-ploc[1,i]
           bestycoord[i]<-ploc[2,i]
         }
         if((pname[i]-1)*20+19 == period) {
           ploc[1,i]<-bestxcoord[i]
           ploc[2,i]<-bestycoord[i]
           toptotals[i]<-0
           indicator[i]<-1
         } else {
           ploc[1,i]<-rnorm(1,mloc[1],1)
           ploc[2,i]<-rnorm(1,mloc[2],1)
         } 
       }
     }
   }
  }

  # 27.Pragmatist move to weighted (alpha) mean of all parties location and their own voter
  pragmatists <- which(ptype==27)
  lpragmatists<-length(pragmatists)
  if(lpragmatists>0) {
   ploc[1,pragmatists]<- rnorm(lpragmatists,mloc[1],0.2)*0.2+pmloc[1,pragmatists]*0.8
   ploc[2,pragmatists]<- rnorm(lpragmatists,mloc[2],0.2)*0.2+pmloc[2,pragmatists]*0.8
  }

  # 28.Raptors
  raptors<-which(ptype==28)
  if (length(raptors)>0) {
        if (! exists("speed") || length(speed) != pn) speed = rep(0,pn)
        close=rep(F,pn)
        if (any(round(speed,3)==0.178)) {
                for (y in raptors) {
                        myspeed=speed
                        myspeed[y]=0
                        for (s in which(round(myspeed,3) == 0.178)) {
                          close[y] = close[y] || (sqrt(sum((ploc[,y] - ploc[,s])^2))< 0.3 &&
                                     ploc[1,y] <= ploc[1,s])
                        }
                }
        }
        angle[raptors]<-(angle[raptors]+(ototals[raptors]>=ptotals[raptors] | close[raptors]) * 
                        runif(length(raptors),pi/2,3*pi/2))%%(2*pi)
        ploc[1,raptors]<-ploc[1,raptors]+0.178*cos(angle[raptors])
        ploc[2,raptors]<-ploc[2,raptors]+0.178*sin(angle[raptors])

        speed = sqrt(colSums((ploc-oploc)^2))   # preserved till next time
  }

        # 29.ShuffleKAN
        ## draw party type 1=sticker, 2=aggregator, 3=hunter, 4=predator, 5=Shufflekan
        TypeNumShuffleKAN <- 29
        ThresKAN <- 0.080*vn
        ThresLargeKAN <- 0.115*vn
        ShuffleKAN <-  which( (ptype==TypeNumShuffleKAN )*(ptotals > ThresLargeKAN) == 1)
        # if the number of supporters is greater than ThresLargeKAN, then this party behaves like 
        # aggregator (but ploc[,ShuffleKAN]<-0.8*pmloc[,ShuffleKAN]) with 80% chance, 
        # or it sticks. 
        if(length(ShuffleKAN)>0) {
                tempKAN <-  (runif(length(ShuffleKAN),0,1) > 0.2 )
                ploc[,ShuffleKAN]<-0.8*pmloc[,ShuffleKAN]* tempKAN - ploc[,ShuffleKAN]*(tempKAN-1)
        }
        # if the number of supporters is less than ThresLargeKAN, then this party 
        # behaves like hunter (but stepsize is 90%) 
        ShuffleKAN<-  which( (ptype==TypeNumShuffleKAN )*(ptotals < ThresLargeKAN) == 1)
        if(length(ShuffleKAN)>0) {
                angle[ShuffleKAN]<-(angle[ShuffleKAN]+(ototals[ShuffleKAN]>=ptotals[ShuffleKAN])*
                        runif(length(ShuffleKAN),pi/2,pi))%%(2*pi)
                ploc[1,ShuffleKAN]<-ploc[1,ShuffleKAN]+ 0.9*pstepsize*cos(angle[ShuffleKAN])
                ploc[2,ShuffleKAN]<-ploc[2,ShuffleKAN]+ 0.9*pstepsize*sin(angle[ShuffleKAN])
        }     
        # if the votes are less than ThresKAN, then LET'S SHUFFLE!
        # JUMP and ESCAPE to a quadrand in which parties have the highest votes per party on average. 
        ShuffleKAN<-  which( (ptype==TypeNumShuffleKAN )*(ptotals < ThresKAN) == 1)
        if(length(ShuffleKAN)>0) {
                ploc[,ShuffleKAN] <-  abs( rbind(rnorm(length(ShuffleKAN)),rnorm(length(ShuffleKAN))) )
                # This will be shifted to the center of gravity, so the party can do this without 
                # knowing the absolute origin (0,0).
                CenterKAN <- mloc
                QuadrandKAN <- matrix(0,4)
                QuadrandKAN[1] <- sum((ploc[1,] >CenterKAN[1])*(ploc[2,] >CenterKAN[2])*ptotals)/
                                        sum((ploc[1,] >CenterKAN[1])*(ploc[2,] >CenterKAN[2]))
                QuadrandKAN[2] <- sum((ploc[1,] <CenterKAN[1])*(ploc[2,] >CenterKAN[2])*ptotals)/
                                        sum((ploc[1,] <CenterKAN[1])*(ploc[2,] >CenterKAN[2]))
                QuadrandKAN[3] <- sum((ploc[1,] <CenterKAN[1])*(ploc[2,] <CenterKAN[2])*ptotals)/
                                        sum((ploc[1,] <CenterKAN[1])*(ploc[2,] <CenterKAN[2]))
                QuadrandKAN[4] <- sum((ploc[1,] >CenterKAN[1])*(ploc[2,] <CenterKAN[2])*ptotals)/
                                        sum((ploc[1,] >CenterKAN[1])*(ploc[2,] <CenterKAN[2]))
                TargetQuadKAN <- which.max(QuadrandKAN)
                if( TargetQuadKAN > 2){
                        ploc[2,ShuffleKAN] <-   -1*ploc[2,ShuffleKAN] 
                }
                if( TargetQuadKAN > 1 & TargetQuadKAN < 4){
                        ploc[1,ShuffleKAN] <-   -1*ploc[1,ShuffleKAN]
                }
                ploc[,ShuffleKAN] <- ploc[,ShuffleKAN]  + CenterKAN
        }



  #
  # store information
  #


  # store previous election totals and locations
  ototals<-ptotals
  umlocch<-(rowMeans(ploc)-rowMeans(oploc))   # unweighted mean change in party location
  oploc<-ploc


 }

}
